// the compiled Rs-DeltaS relation
// created: March 15, 2021 by Yazhou SUN
// lasted modified: May 06, 2022

double Rs(double rs, bool opt);
double DS(double ds, bool opt);

struct tRs{
  string name;
  double rs, drs;
  double ds;
  tRs(const string &n, double rs_, double drsm, double drsp, double ds_, bool opt1 = true, bool opt2 = true){
    name = n;
    if(fabs(rs_-(drsm+drsp)/2.) > 0.3){
      cout << name << ": large discrepancy happened." << endl;
      cout << "rs_: " << rs_ << " drsm: " << drsm << " drsp: " << drsp << endl;
      getchar();
    }

    rs = Rs(rs_, opt1);
    drs = fabs(Rs(drsm, opt1)-Rs(drsp, opt1))/2.;
    ds = DS(ds_, opt2);
  } // end the constructor
  void Print(){
    cout << name << " rs: " << rs << "(" << drs << ")," << " ds: " << ds << endl;
  }
};

void rsdsCmpHFCp1(){
  vector<tRs> vn, vp, ve;

  vn.push_back(tRs("Be10-n", 0.35, -0.7, 1.4, 11.1));
  vn.push_back(tRs("Ar46-n", 4.5, 1.4, 7.7, 12.55));
  vn.push_back(tRs("Li9-n", 4.75, 3.75, 5.85, 13.02));
  vn.push_back(tRs("Si38-n", 4.8, 4.0, 5.6, 9.27));
  vn.push_back(tRs("Si40-n", 5.0, 4.0, 6.1, 7.81));
  vn.push_back(tRs("Si36-n", 3.0, 2.45, 3.5, 10.8));
  vn.push_back(tRs("C15-n", 7.45, 6.4, 8.5, 6.51));
  vn.push_back(tRs("C19-n", 7.9, 5.8, 10.0, 2.45));
  vn.push_back(tRs("C20-n", 8.2, 6.4, 10.0, 2.80));
  vn.push_back(tRs("Ni57-n", -3.8, -4.3, -3.3, 11.75, 0,0));
  vn.push_back(tRs("C12-n", -3.55, -4.05, -3.0, 13.74, 0,0));
  vn.push_back(tRs("O16-n", -1.75, -2.5, -0.95, 10.75, 0,0));
  vn.push_back(tRs("Mg22-n", -7.2, -7.9, -6.4, 6.66, 0,0));
  vn.push_back(tRs("Ar34-n", -5.65, -7.95, -3.8, 6.62, 0,0));
  vn.push_back(tRs("Ca36-n", -10.0, -10.55, -9.4, 4.8, 0,0));
  vn.push_back(tRs("Ar32-n", -9.75, -10.55, -8.9, 3.16, 0,0));
  vn.push_back(tRs("S28-n", -8.2, -7.4, -8.9, 3.26, 0,0));
  vn.push_back(tRs("Si24-n", -6.15, -7.2, -5.1, 4.1, 0,0));
  vn.push_back(tRs("C10-n", -3.6, -4.05, -3.0, 4.4, 0,0));

  vp.push_back(tRs("O14-p", 2.2, 5.1, -0.7, 7.38));
  vp.push_back(tRs("Si24-p", 3.2, 2.2, 4.3, 7.86));
  vp.push_back(tRs("Ca36-p", 3.75, 2.7, 4.8, 8.5));
  vp.push_back(tRs("C9-p", 8.4, 7.1, 9.7, 10.8));
  vp.push_back(tRs("C9-p", 3.75, 2.2, 5.3, 11.05));
  vp.push_back(tRs("Be8-p", 4.8, 3.0, 6.65, 11.15));
  vp.push_back(tRs("S28-p", 6.4, 4.6, 8.2, 7.06));
  vp.push_back(tRs("C12-p", -3.7, -4.25, -3.2, 17.35, 1,0));
  vp.push_back(tRs("O16-p", 1.4, 0.4, 2.4, 15.25, 0,0));
  vp.push_back(tRs("Be10-p", -5.35, -5.9, -4.8, 7.37, 0,0));
  vp.push_back(tRs("Mg28-p", 0.6, 0.1, 1.1, 10.75, 0,0));
  vp.push_back(tRs("Si36-p", -6.15, -6.65, -5.6, 7.00, 0,0));
  vp.push_back(tRs("Si38-p", -8.45, -8.9, -7.9, 5.53, 0,0));
  vp.push_back(tRs("Si40-p", -8.45, -10.25, -6.65, 4.04, 0,0));

  ve.push_back(tRs("Ca40-e", -0.65, -1.95, 0.65, 14.76));
  ve.push_back(tRs("P31-e", 0.15, -1.1, 1.15, 18.95, 1,0));
  ve.push_back(tRs("O16-e", -1.1, -2.7, 0.9, 18.0, 1,0));
  ve.push_back(tRs("C12-e", -1.9, -0.9, -3.0, 17.5, 1,0));
  ve.push_back(tRs("Zr90-e", -3.0, -4.25, -1.7, 18.05, 1,0));
  ve.push_back(tRs("V51-e", -4.8, -5.8, -3.75, 17.7, 1,0));
  ve.push_back(tRs("Li7-e", 0.1, 1.6, -1.5, 13.9, 0,0));
  ve.push_back(tRs("Si30-e", -1.2, -2.5, 0.1, 13.8, 0,0));
  ve.push_back(tRs("Ca48-e", -2.25, -3.3, -1.2, 11.85, 0,0));
  ve.push_back(tRs("Pb208-e", -3.55, -4.8, -2.2, 15.35, 0,0));

  const int nn = vn.size(), np = vp.size(), ne = ve.size();
  double rs_n[nn], ers_n[nn], ds_n[nn], eds_n[nn];
  double rs_p[np], ers_p[np], ds_p[np], eds_p[np];
  double rs_e[ne], ers_e[ne], ds_e[ne], eds_e[ne];
  for(int i = 0; i < nn; i++){
    rs_n[i] = vn[i].rs; ers_n[i] = vn[i].drs; ds_n[i] = vn[i].ds; eds_n[i] = 0.;
    // vn[i].Print();
  }
  for(int i = 0; i < np; i++){
    rs_p[i] = vp[i].rs; ers_p[i] = vp[i].drs; ds_p[i] = vp[i].ds; eds_p[i] = 0.;
    // vp[i].Print();
  }
  for(int i = 0; i < ne; i++){
    rs_e[i] = ve[i].rs; ers_e[i] = ve[i].drs; ds_e[i] = ve[i].ds; eds_e[i] = 0.;
    // ve[i].Print();
  }
  TGraphErrors *gn = new TGraphErrors(nn, ds_n, rs_n, eds_n, ers_n);
  TGraphErrors *gp = new TGraphErrors(np, ds_p, rs_p, eds_p, ers_p);
  TGraphErrors *ge = new TGraphErrors(ne, ds_e, rs_e, eds_e, ers_e);
  gn->SetMinimum(0.15); gn->SetMaximum(1.3);
  gn->SetTitle("");
  TAxis *ax = gn->GetXaxis(), *ay = gn->GetYaxis();
  ax->SetRangeUser(-30, 24);
  ax->SetTitle("\\DeltaS (MeV)"); ay->SetTitle("R_{s}=\\sigma_{exp}/\\sigma_{th}");
  ax->SetTitleOffset(1.0);

  ax->SetLabelSize(0.05); ax->SetTitleSize(0.05);
  ay->SetLabelSize(0.05); ay->SetTitleSize(0.05);
  ax->SetTitleOffset(0.85); ay->SetTitleOffset(0.95);
  ax->CenterTitle(); ay->CenterTitle();

  double markersize = 2, linewidth = 1, col = 17;
  gn->SetMarkerStyle(20); gn->SetMarkerSize(markersize);
  gn->SetMarkerColor(21); gn->SetLineColor(21); gn->SetLineWidth(linewidth);
  gp->SetMarkerStyle(22); gp->SetMarkerSize(markersize);
  gp->SetMarkerColor(40); gp->SetLineColor(col); gp->SetLineWidth(linewidth);
  ge->SetMarkerStyle(21); ge->SetMarkerSize(markersize);
  ge->SetMarkerColor(45); ge->SetLineColor(col);  ge->SetLineWidth(linewidth);
  gStyle->SetEndErrorSize(9);
  TLegend *lg1 = new TLegend(0.105,0.105,0.45,0.3);
  lg1->AddEntry(ge, "(e,e' p): \\DeltaS=S_{p}-S_{n}", "p");
  lg1->AddEntry(gn, "n-removal: \\DeltaS=S_{n}-S_{p}", "p");
  lg1->AddEntry(gp, "p-removal: \\DeltaS=S_{p}-S_{n}", "p");
  lg1->SetBorderSize(0);

// systematics of gade-tostevin-2014
//-------------------------------------------------------


  const int n = 2;
  const double xh[n] = {-30., 30.}, yh[n] = {1.23, 0.4};
  TGraph *gHatch = new TGraph(2, xh, yh);
  gHatch->SetLineColor(kWhite);
  gHatch->SetLineWidth(-3200);
  gHatch->SetFillStyle(3004);
  gHatch->SetFillColor(1);
  gHatch->SetMinimum(0.15); gHatch->SetMaximum(1.3);
  gHatch->SetTitle("");
  // central and envelope lines //
  const double dy = 0.18;
  const double xu[n] = {-30., 30.}, yu[n] = {1.23, 0.4};
  const double xc[n] = {-30., 30.}, yc[n] = {1.23-dy, 0.4-dy};
  const double xd[n] = {-30., 30.}, yd[n] = {1.23-2.*dy, 0.4-2.*dy};
  TGraph *gc = new TGraph(2, xc, yc); gc->SetLineStyle(1);
  TGraph *gu = new TGraph(2, xu, yu); gu->SetLineStyle(2);
  TGraph *gd = new TGraph(2, xd, yd); gd->SetLineStyle(2);
  // TAxis *ax = gHatch->GetXaxis(), *ay = gHatch->GetYaxis();
  // ax->SetRangeUser(-30, 24);
  // ax->SetTitle("\\DeltaS (MeV)"); ay->SetTitle("R_{s}=\\sigma_{exp}/\\sigma_{th}");
  // ax->SetTitleOffset(1.0);
  //
  // ax->SetLabelSize(0.05); ax->SetTitleSize(0.05);
  // ay->SetLabelSize(0.05); ay->SetTitleSize(0.05);
  // ax->SetTitleOffset(0.85); ay->SetTitleOffset(0.95);
  // ax->CenterTitle(); ay->CenterTitle();
  // gStyle->SetEndErrorSize(9);

  // r0 = HF //
  TGraphErrors *gpHF = new TGraphErrors();
  gpHF->SetPoint(0, -1.796642, 0.488437); gpHF->SetPointError(0, 0., 0.042167); // C12
  gpHF->SetPoint(1, 13.314627, 0.418371); gpHF->SetPointError(1, 0., 0.063210); // C13
  gpHF->SetPoint(2, 12.870917, 0.424864); gpHF->SetPointError(2, 0., 0.028120); // C14
  gpHF->SetPoint(3, 20.136266, 0.431190); gpHF->SetPointError(3, 0., 0.042345); // C15
  gpHF->SetPoint(4, 18.302800, 0.290464); gpHF->SetPointError(4, 0., 0.031742); // C16
  gpHF->SetMarkerStyle(20); gpHF->SetMarkerSize(2);
  gpHF->SetMarkerColor(kRed); gpHF->SetLineColor(kRed);
  gpHF->SetLineWidth(3);

  TLegend *lg = new TLegend(0.505,0.65,0.895,0.86);
  lg->AddEntry(gpHF, "R_{s} for ^{x}C-p from ETF", "p");
  lg->SetBorderSize(0);

  TCanvas *c = new TCanvas("c", "Rs-DeltaS", 800, 600);

  gn->Draw("Ap");
  ge->Draw("samep");
  gp->Draw("samep");
  gHatch->Draw("sameC");
  lg1->Draw();
  gc->Draw("lsame"); gu->Draw("lsame"); gd->Draw("lsame");
  TF1 *fpol1p = new TF1("fpol1p", "pol1", -30., 25.);
  fpol1p->SetLineColor(gpHF->GetLineColor());
  // gpHF->Fit("fpol1p");
  gpHF->Draw("samep");
  lg->Draw();

  TPaveText *pt = new TPaveText(-2.739644,0.5468334,-0.3535978,0.609409,"br");
  pt->SetBorderSize(0);
  pt->SetFillColor(0);
  pt->SetTextColor(2);
  pt->SetTextFont(42);
  TText *pt_LaTex = pt->AddText("^{12}C");
  pt->Draw();

  pt = new TPaveText(14.43429,0.4554687,16.74857,0.5153646,"br");
  pt->SetBorderSize(0);
  pt->SetFillColor(0);
  pt->SetTextColor(2);
  pt->SetTextFont(42);
  pt_LaTex = pt->AddText("^{13}C");
  pt->Draw();

  pt = new TPaveText(9.205714,0.4305121,11.60571,0.4929036,"br");
  pt->SetBorderSize(0);
  pt->SetFillColor(0);
  pt->SetTextColor(2);
  pt->SetTextFont(42);
  pt_LaTex = pt->AddText("^{14}C");
  pt->Draw();

  pt = new TPaveText(21.46286,0.4517253,23.77714,0.5141168,"br");
  pt->SetBorderSize(0);
  pt->SetFillColor(0);
  pt->SetTextColor(2);
  pt->SetTextFont(42);
  pt_LaTex = pt->AddText("^{15}C");
  pt->Draw();

  pt = new TPaveText(19.40571,0.2158854,21.72,0.2782769,"br");
  pt->SetBorderSize(0);
  pt->SetFillColor(0);
  pt->SetTextColor(2);
  pt->SetTextFont(42);
  pt_LaTex = pt->AddText("^{16}C");
  pt->Draw();

  c->Print("rsCmpHFCp1.pdf");
  // c->Print("/home/asia/Lab/cnok/cnok/docs/rs22/rsCmpHFC.pdf");
  // c->Print("rsCmp.eps");
} // end the main function

// -----------
double Rs1(double rs){
  static const double kRs = 0.8/(11.+9.82), bRs = 0.76-2.2*kRs;
  return kRs*rs+bRs;
}
double Rs2(double rs){
  static const double kRs = 0.8/(11.+9.82), bRs = 0.5+3.85*kRs;
  return kRs*rs+bRs;
}
double DS1(double ds){
  static const double kds = 16./10.5, bds = -10.-13.02*kds;
  return kds*ds+bds;
}
double DS2(double ds){
  static const double kds = -16./10.5, bds = -3.-17.7*kds;
  return kds*ds+bds;
}
double Rs(double rs, bool opt){
  if(opt) return Rs1(rs);
  return Rs2(rs);
}
double DS(double ds, bool opt){
  if(opt) return DS1(ds);
  return DS2(ds);
}
